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1. Introduction 

Solving time dependent or viscoelastic problems for a 
homogeneous isotropic material can be involved and tedious. 
Extending this to nonhomogeneous and anisotropic materials 
such as layered fiber reinforced composite materials can be 
nearly impossible for closed form solutions. However, with 
numerical methods the designer or engineer of these 
materials can predict, with reasonable accuracy, the 
viscoelastic response without doing actual creep tests on 
each possible laminate. 

The overall criteria for an acceptable viscoelastic 
numerical method is one that will be stable for large time 
steps, converge to the correct answer, and not take a 
tremendous amount of computer resources and time. In 
addition to these, the program that will use the method, 
will have to run on a microcomputer which further restricts 
the maximum run time and total memory. This will allow easy 
access to the program for design engineers and will make the 
design process, with its many ’what if’ conditions and 
numerous rerunning proceed faster. 

This report will examine various numerical methods that 
have been used in solving numerical viscoelastic methods. A 
new method, called the Nonlinear Differential Equation 
Method CNDEMD , which is based on the prony series, will be 
introduced and compared with the current methods. The later 
part of this report will deal with the actual implementation 
and verification of the NDEM method. 


PRECEDING PAGE BLANK NOT FILMED 
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2. Previous Work at VPI 

The concept of predicting the viscoelastic response in 
any general laminate has been previously investigated by 
others at Virginia Polytechnic Institute and State 
University. Dillard, et al,Cl,2] first formally proposed 
using known unidirectional material properties C obtained 
experimentally? of a composite lamina to predict the 
nonlinear viscoelastic response of any general laminate 
constructed from the same material by numerical methods. 

They examined the graphite/epoxy T300/934 composite system 
and closely predicted the response of various general 
laminate composites. Others, Tuttle [33 and Heil [43, have 
also used this basic concept to closely predict the response 
of other graphite/epoxy systems. 

The numerical solution method used by Dillard C23 was 
based on classical lamination theory, with time incremented 
in a step fashion. The solution scheme first calculates the 
static stress and then begins the time step increments. The 
strain state is determined at t+At, using the stress state 
at time t and the viscoelastic constitutive equation for 
that particular material. The stress state is assumed to be 
constant throughout the time step from t to t+At. The new 
ply stresses are then determined at t+At based on the 
current creep strains and the applied mechanical load. This 
cycle is repeated, with the new stresses substituted back 
into the nonlinear compliance functions, until the stresses 
converge. A new time step is then taken and the processes is 
repeated. The algorithm for calculating creep strains is 
similar to the classical lamination theory method of 
calculating the strains due to thermal loads. This 
procedure was implemented on an IBM mainframe computer and 
was called VI SLAP C VIScoelastic LAmi nation Program? by 


< 
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There were three major difficulties with VISLAP and its 
numerical method. One. the basic algorithm of substituting 
old stresses back into the nonlinear compliance functions, 
and repeating the solution process until all stresses 
converge can have stability problems. This algorithm of 
successively substituting an unknown variable into a set of 
equations until convergence, called the Gauss-Seidel or 
successive substitution method, is not unconditionally 
stable. For example, if the coefficient matrix, CC3, in the 
following set of equations, represented in matrix form, 

[C3<x> = [B] Cl. ID 

is not positive definite then it will not converge [53. In 
some laminate cases, predominantly two fiber angle 
laminates, VISLAP will be unstable for this reason. 

The second difficulty with VISLAP concerns the large 
time step size necessary to reach a solution of problems 
covering long time spans. If the time step is sufficiently 
large, stability problems will arise. VISLAP basically uses 
a first order forward integrating method, called the Euler 
Method [63, to solve for the creep strains at each step, 
which will have a maximum step size to remain stable. 

In conjunction with the time step size problem is the 
third difficulty with VISLAP, the actual computer time and 
computer memory space needed for a solution grows 
exponentially with each additional time step. As each time 
step is taken, the creep strain must be recalculated over 
the entire time span back to the initial start time. This 
requires that all stresses at each time step must be stored 
and used for calculations at the next time step. This 
recalculation of the creep strain integral at each time step 
becomes more time consuming with each additional step. In 
order to minimize the computer solution time and memory, 
step sizes are increased in a logarithmic manner as the 
solution progresses. However, as stated above, large time 
steps can cause stability problems. 
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In order to overcome some of the problems in VI SLAP but 
still retain its ability to calculate the complex, time 
dependent stress and strain state of an orthotropic 
composite laminate, various common numerical solution 
techniques will be investigated in the following sections. 

In addition to the those methods, a new method will be 
presented which resolves all the problems dealing with 
stability and solution time length. 
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3. Direct Iteration of the Volterra Integral 


Viscoelastic problems naturally fall into the broad 
class of mathematical problems called convolution integral 
equations, of which the Volterra integral equation of the 
second kind is the most common. The general form of the 
Volterra equation is 

X 

uCxD = fCxD + Xj kCx, tD uCtD dt C3. ID 

a 

where uCxD is the unknown function and fCxD, kCx.tD, and X 
are known functions or constants. By simply changing the 
variable and function names and forms, the well known 
hereditary integral in viscoelasticity [7] becomes evident. 

t 

«CtD = oCoDDCtD + J DCt-iO dr C3.2D 

o 

where sCtl is the total strain, DCtD is the compliance 
function* and crCrD is the stress function. This form is for 
a single homogeneous material. For a material made from 
multiple homogeneous layers, i . e. , composite laminates, the 
total strain rCtD will be dependent on each of the stress in 
each of the layers. The strain can be written in terms of 
stress as 

t 

tfCoCtO.O ■ oCoDDCtD + J DCt-iO ' - dr C3. 3D 

o 

which is a Volterra integral of the second kind. 

A simple example of such a system would be a one 
dimensional laminate material that is constructed from two 
parallel materials as illustrated in Fig. 3.1. The two 
materials have different compliance functions and the 
complete laminate is under a constant laod. In this example 
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one material will be elastic and the other viscoelastic with 
a ^ and a^ as the stresses in material 1 and 2, respectively 
and as the total applied stress. If DCO represents any 
time dependent compliance function then 

p l da Ct3 

c CO = o' Co3DC t3 + DCt-r3 dr C3.43 

a a J dr 

o 


By using the relationships 

e Ct3 + e Ct) = eCO 
1 2 

a CO + a C O = a C O = a 
12 0 0 

E *:Ct3 * O' CO 

1 l 

E eCO * O’ CO 

2 2 

it can be shown that 


C 3. 5a3 
C 3. 5b5 
C3. 5c3 
C3. 5d3 


o- Ct3 

l 

E 


i 


a Ct3 
2 

E 


2 


t 

J DC t-O 

o 


da Ct3 
2 

dr 


dr 


+ o- Co3DCt3 
2 


C3. 63 


or 


a Ct3 
2 



r da Ct3 

o' E -EE DCt~r3 1- dr 

O 2 1 2J dr 

o 

- E E o- Co3DCO C 3. 73 

12 2 


This can be further simplified by integrating by parts to 
give 


O' CO = 

2 


a E 

0 1 

E +E 

1 2 


E E 
1 2 

E +E 
1 2 


V 

I 


dDCt-O 

dr 


a Ct3 
2 


dr 


C 3. 83 


This form can be more easily evaluated since DCt3 is usually 
given and its derivative can be calculated directly, where 

a Ct3 is not known and its derivative is difficult to find. 

2 

Equation 3.8 is in the standard convolution Vol terra 
integral form which has been studied in detail by others 
[8-103 from a mathematical point of view. Once a ^ CO is 


Page 9 


known, the total strain cCtJ can easily be calculated from 
Eq. 3.5. It should be noted that the one dimensional 
example presented a very simplified case and for a more 
natural multidimensional material the equation would not 
only be more complex, but there would be several coupled 
equations and not just one. However, to understand the 
basic principles and difficulties in solving the Volterra by 
numerical methods, the given example will be examined. 

A closed form solution of Eq. 3.8 is possible for 
certain compliance functions, DCt3, such as a linear dashpot 
model, DCt3 = t/^i, or a Kelvin element model, DCt3 ■ 

1 -exp[ -tCE/juD ] . However, compliance functions with 
solutions are scarce and are found for only simple 
functions. One important function that is widely used in 
linear viscoelastic analysis, and does not have a closed 
form solution, is the power law equation 

DC t3 = mt n C3. 93 

where m and n are constants. Since closed form solutions 
are difficult to obtain and limited to certain compliance 
functions, numerical methods need to be applied to obtain 
most solutions. 

Four concepts to be considered when employing numerical 

methods are convergence, error, stability, and solution 

time. The solution time length becomes especially critical 

when dealing with convolution Volterra integrals, due in 

part because dCt3 and dDCt-iO/dr continually change with 

each new time step. This requires the complete integral to 

be recalculated for each new time step. Unlike standard 

integrals, past results can not be used to calculate future 

points, but the total integral, from t C t =03 to the current 

o 

time, t, must be recalculated. At long times, i.e. large 
number of time steps, this method can require a tremendous 
amount of time and computer memory storage. If, however, 
the time steps can be varied, such as short steps at the 
start where the function is changing rapidly and long steps 



Page 1 O 


towards the end where the function is changing slowly, then 
this method can be economical. 

Convergence is generally not a problem f or a non-singular 
kernel or compliance function, DCtD. However it does become 
a concern if the kernel is not well behaved or is singular. 

The power law C equation 3. 9D , which is used extensively in 
viscoelasticity, is classified as weakly singular, meaning 
the derivative at some point is singular or undefined Cat 
zero for the power lawD. The solution of the integral can 
converge with weakly singular functions if the time steps 
around the weakly singular point are sufficiently small. 
Convergence of the power law and its associated problems 
will be demonstrated with an example later in this section. 

Stability or numerical oscillations can occur in the 
solution of numerical problems. Even if the problem seems 
to converge and the error is small , it could diverge after a 
certain time step or step size. Two common causes of 
stability problems are: ID the numerical precision of the 

computer or code, which leads to round off errors and 
truncation, and 2D the time step size. Generally the 
precision of the computer is not a problem or can be solved 
by upgrading to a better computer or programming language. 

On the other hand, most numerical solution techniques have a 
limit on the time step size before stability becomes a 
concern. All forward or explicit numerical integration 
techniques, which are generally used for the convolution 
integral, are not absolutely stable for all time step sizes 
[6,8,9,101. This is a serious concern with viscoelastic 
analysis since increasing time steps are necessary to reduce 
the computer calculation time and memory size, as explained 
in the preceding paragraphs. 

Error is associated with the accuracy of the computer 
and the algorithm used to solve the problem. Various 
algorithms have been developed for the solution of 
convolution integral equations which include, in ascending 
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order of accuracy, Euler, Modified Euler or trapezoidal, 
Simpson rule with trapezoidal end, Simpson rule with ^ rule, 
and Runge-Kutta. The higher order methods take more time 
for each time step but the accuracy is generally higher and 
larger time step sizes are possible. The trapezoidal 
algorithm will be presented in detail to demonstrate how the 
Volt err a Integral can be solved numerically. Other methods 
are similar and will not be presented. However, the 
solution of the example problem presented earlier by all 
methods mentioned above will be compared at the end of this 
section. 

An approximation for the convolution integral can be 
written 


J KC^-tO 

o 


oCtD dr = h \ w KCt.-TO oC r.) 

j*o 1 J J 


h 



w. . K oC T .5 
>-J J 


i= 0, 1 , 2, . . N 


C3. io:> 


where h is the step size, K. . * dDCt.-rD/dr, and w are 

IJ v j ' Vj 

the weights for the appropriate integration rule. For 
example, the weights for the trapezoidal method are w. = 

LO 

w..= and w ™ 1. All or “the preceding weights assume 
equal step sizes. In this manner the first few steps of Eq. 
3.8 for the trapezoidal method are 


a E 

act =03 = ? p 

2 o E + E 


1 2 
ct E E E h 


V C. Ci & II »» 

Ct) = -= ° w K O' Ct 3 + w K O' Ct D 

i E + E E + E I 10 10 2 o ii 11 2i 

1 2 1 2 "» J 

ct E E E h 

= + - ~Tv - K O' Ct D + - K o' Ct D I 

E + E E + E I 2 14 2 o 2 11 2 1 

1 2 1 2 t J 
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aCi) 

2 2 


o' Ct ) 

2 V 


o* E E E h r 

Ol 12 | 1 

S + E E + E I 2 

12 1 2 *• 


K oCt ) + K c^CtD 

20 2 O 2 12 1 


+ - E oCt ) 
2 22 2 2 


] 


O' E 


E E 


V C, E. C. . I* -I 

•°— | r- + — - 1 * — 5 I K «7(t ) + K ffCt)] 

+ E E + E 2 I to 2 o u z 1 

1 2 1 2 L J 

i v 1 

+ > K o- Ct D C3. lla-dD 
. L t j 2 j 


j = 1 


In each of these steps the unknown stress, o’ Ct D, can be 

2 v 

factored out and solved for by manipulating the equation 
algebraically. However, if the kernel KCtD is nonlinear in 
terms of stress, then all i nonlinear equations would need to 
be solved simultaneously. This quickly becomes prohibitive 
since there will be thousands of time steps in a typical 
problem, which translates to solving thousands of nonlinear 
equations simultaneously. Similar relationships to equation 
3.11 can be constructed for other integrating schemes. For 
higher order methods such as the Simpson rule or Runge — 
Kutta, a starting procedure needs to be used which should be 
of the same order of magnitude in accuracy. Various 
starting techniques can be found in the literature [11-143. 

To evaluate the use of the Vol terra integral for 
viscoelastic materials the one dimensional example described 
at the beginning of this section CFig 3. ID will be used. 

Two different but common compliance functions, DCtD, were 
chosen to be examined, a dashpot, DCtD = t//u, where ^ is the 
viscosity constant of the dashpot and a power law, DCtD = 
mt n , where m and n are assumed given. 

The dashpot function has a exact solution to Eq. 3. 8, 
which will be used to verify the numerical results. 


cr 

2 


O' 

o 


E 

2 


E + E 


1 2 


-\t 


e 


C3. 12D 


where X = E E /uCE + E D and a is constant. Five different 

12' 1 2 O 
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integrating techniques were used to solve the example 
problem: Euler, trapezoidal, trapezoidal with 3-point 

starting technique, Simpson with 3-point starting technique 
and trapezoidal rule for even last data point, and Simpson 
with 3-point starting technique and 3/8 rule for the even 
last data point. The results are shown on Fig 3.2. Even 
though the time step, h, was large, all but the Euler method 
are within acceptable accuracy limits. 

The second compliance function to be examined, the 
power law, has no closed form solution to compare with the 
numerical results. However, by examining the results of 
various integrating techniques the solution can be deduced. 
The same five integrating techniques used for the dashpot 
test case were also used for the power law C constants m = 5 
and n = 0. 2D and the results are shown in Fig. 3.3. The 
time step was h = 0.1, two magnitudes smaller then for the 
dashpot example, but unlike the dashpot results the power 
law results vary and even oscillate. If the step size is 
reduced, the solution tends to converge to smaller values 
CFig. 3. 4D and it becomes evident that the time step size 
affects the solution convergence. The solution does seem to 
slowly approach a limiting value as h ** 0.0. 

The solution of the power law function is inaccurate 
because it is a weakly singular function at zero. The 
derivative of the power law at zero is infinity and the 
derivative changes rapidly for small values of time. This 
requires very small time steps, C= 10 near the origin for 
any of the numerical integration techniques to converge. 
However, with small time steps, the time required to solve 
the problem increases tremendously which then limits the 
time span. The time step size can be increased as the time 
becomes larger but there will be a upper limit on step size 
before stability difficulties develop. 

In conclusion, the direct numerical integration of the 
Vol terra integral for linear viscoelastic problems is not 
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recommended. The biggest difficulty was the long run times 
necessary for any numerical solution to converge when using 
the power law compliance function. This was caused by the 
weakly singular nature of the power law. Other difficulties 
would be the inclusion of nonlinear stress effects, thus 
creating a large number Con the order of hundreds} nonlinear 
equation that would need to be solved simultaneously. It 
should also be noted that the above difficulties would be 
magnified for multidimensional materials such as orthotropic 
composite materials. 
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4 . Prony Series in Modeling Linear Viscoelastic Response 


The Prony series is a method to model viscoelastic 
response that is derived from a series of Kelvin elements. 
This series can best be understood by first looking at a 
single Kelvin element, which has a spring and dashpot in 
parallel as shown in Fig. 4.1. The Kelvin element needs two 
parameters to describe its response to a given load or 
displacement, the spring constant, E .and the dashpot 
viscosity, The load or stress, a ^ and the strain, £, can 

be related by summing the stress in both the spring and 
dashpot. 

a + a , = a C 4 . 1 D 

s d s 

Substituting the constitutive equations for a spring and 
dashpot gives 

cE + cu = a C 4 . 2D 

^ o 

Solving for £, and assuming is constant will give 

This can be generalized with a series of Kelvin elements as 


£ 




C 4. 4D 


where n is the total number of Kelvin elements in the 
series. As a further generalization, a single spring can be 
placed in series with the Kelvin elements such that 



C4. 5D 


where is the spring constant in the single spring. 
Equation 4.5 is referred to as a Prony series. Prony series 
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can model the creep of a viscoelastic material accurately if 

the retardation times Cr = p/E} of the individual Kelvin ^ 

elements are properly spaced. Since one Kelvin element 

influences the strain over about 1 1/2 decades of time, the 

retardation times should be spaced about one per decade of 

time that is being modeled. * 

One advantage of the Prony series is its ability to 
accurately represent any data over any time span if enough 
elements are used. This is especially useful if that data 
is not uniform or does not conform to any general curve < 

shape. Of course, this is also a disadvantage since a large 
number of material properties, two for every element, are 
required. In a contrast, the linear power only has 3 

parameters, e — e + mt n » to describe the viscoelastic ( 

strain response. Another major difference between the Prony 

series and the power law is the extrapolation of creep 

response outside the actual collected data range. The Prony 

series is derived or fitted only to actual data and after 1 

the last data point the series stops. The power law is also 

derived from actual data but after the last data point the 

equation still indicates or predicts a change in creep over 

time. Although, prudent engineering prohibits the use or 

extrapolation of results past actual collected data, it is 

still useful to understand the expected creep response or 

trend. 

One of the most important advantages of the Prony 
series is that each Kelvin element can be solved 
independently as a differential equation Csee equation 4. 2Z> 
and then the solutions can be summed together. A 
differential equation in the form of equation 4.2, allows 
the use of common and well understood numerical methods for 
solving differential equations. Since the problem has been 
transformed to solving differential equations and not a 
convolution integral Ci.e. , the Vol terra Integral} the 
solution techniques are simpler and easier to implement on a 
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computer. The results of each time step no longer need to 
be stored and reused to calculate future creep steps like 
the convolution method requires. All information needed to 
take another time step is available in the current solution 
of each differential equation. 

The concept of using Kelvin elements and their 
respective differential equations to solve viscoelastic 
problems was presented by Zienkiewicz, et al, [15,16]. They 
used the differential equation formulation in conjunction 
with finite elements method to successfully solve 
geometrically complex problems. The constant stress 
solution, Eq. 4. 3, was used to develop a solution technique. 
By taking a small time step, Eq. 4.3 can be written as 



where ). is the strain from the previous time step 
solution and At is the current time step size. If the 
stress is constant for all time steps, Eq. 4.6 will give an 
exact answer to Eq. 4.2. However in most practical problems 
the stress is constantly changing due to relaxation, 
temperature changes, load changes, etc. If the time step is 
small and the stress changes gradually, then Eq. 4.6 gives 
accurate results as shown by Zienkiewicz. 

In order to describe the viscoelastic response over 
long periods of time, Kelvin elements with different 
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relaxation times may be combined in series as shown in Fig. 
4.2. Each Kelvin element is described by a differential 
equation and the solution of which can be written in the 
form of Eq. 4.6. These solutions can then be summed 
t oget her to gi ve 


e 

T 



+ 




C 4. 7D 


( 


where l is the total number of Kelvin elements in the 
series. 

Solution techniques based on Eq. 4.7 have been widely 
used for stress analysis of linear isotropic materials for 
limited time spans [16-18]. There are three main 
deficiencies with Eq. 4.7 formulation. First, only linear 
viscoelastic materials can be analyzed, whereas many of 
today’s materials, specifically plastics, are nonlinear. A 
nonlinear viscoelastic material will have a different 
compliance and rate of change of compliance at different 
stress levels. Since Eq. 4.7 does not account for these 
nonlinearities, the numerical results will possibly not 
agree with actual experimental results. Some researchers 
[17] have extended Eq. 4.7 to include nonlinear effects with 
limited success. Second, the time step size has an upper 
limit at which the numerical solution technique will become 
unstable since the equation is a forward difference or an 
explicit method. Only an implicit numerical method can be 
’unconditionally stable’ for all time step sizes [63. Limiting 
the time step size in a viscoelastic problem, which can span 
many decades of time, is a concern since large time steps 
become necessary toward the end of the problem. Third, Eq. 

4.7 is only a first order numerical solution technique, 

! 

commonly referred to as the ’Euler Method’, for differential 

i 
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equations. To increase the accuracy and/or decrease the 
number of steps necessary, a higher order solution technique 
should be employed. 

Other general problems with Eq. 4. 7 are the constant 
stress assumption at each time step and the difficulties of 
using it with orthotropic materials. However, the two 
advantages of not having to store all past results and not 
having to recalculate strain at previous time steps for 
every new time step taken overshadows the disadvantages. 

The following section will present a method to extend the 
Prony series method to solve orthotropic, nonlinear 
viscoelastic problems for long time spans. 
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5. Nonlinear Differential Equation Method with the 

Prony Series 

The basic concepts of Kelvin elements and Prony series 
presented in section 4 will be utilized and extended to 
include nonlinear effects, viscoelastic orthotropic 
materials, and unconditionally stable time steps. It was 
these three difficulties that limited the use of the 
Zienkiewicz differential equation method for viscoelastic 
anal ysis. 

The basic differential equation for a single Kelvin 
element can be written as Csee Eq. 4.2 and Fig 4. ID 

^ * -T- ^ - -5T * cs ' i:> 

Where D is the compliance of the spring Cl /ED and X is the 

retardation time Cp/ED. Both the compliance and retardation 

time are considered known and can be obtained from the Prony 

series used to describe the viscoelastic response Csee Eqs. 

4.4 and 4. 5D . A single differential equation of the form 

Eq. 5.1 can be used for each term in a given Prony series. 

Up to this point only one material property has been 

dealt with at a time. However, all materials are defined, 

as a minimum, by at least two material properties which need 

to be considered simultaneously. Isotropic materials are 

generally described by Youngs’ Modulus and Poisson’s ratio, 

whereas orthotropic materials have four independent material 

in-plane properties which are commonly referred to as the 

fiber direction stiffness CE D, the transverse direction 

n 

stiffness CE D, the shear modulus CE or G D , and 

ZZ <50 1Z 

Poisson’s ratio of the fiber direction to transverse 
direction Cv^D. In condensed matrix form, these properties 
relate stress and strain as 
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where s and are strain and stress, respectively, in the 

fiber direction, e and a are strain and stress in the 

2 2 

transverse direction, respectively, and c and <7^ are 
shear strain and stress, respectively. The matrix 
containing E , E , v , and E is referred to as the 
compliance matrix CS] which can be written as 
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where S , S , S , and S are the four independent 

11 12 22 <5*5 r 

properties needed to characterize an orthotropic material. 
These four terms will be referred to as S where q goes from 
1 to 4, such that S=S,S =*S,S = S , and S * S . 

11 1 12 2 22 9 (30 4 

This numbering convention becomes necessary, as will be 
seen later, to differentiate these orthotropic compliance 
terms from the rotated compliance matrix terms, which will 
use the double subscripts CS 

The viscoelastic portion of each of the unrotated, S 
-terms can be described by a Prony series. The general form 
is 


S 

<? 


n r 

1=4 
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1 = 1 ,2,3. . n 
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where ,D is the compliance coefficient variable for the l 
1 <* th 

Kelvin unit in the q direction and is the retardation 


th 


Both ,D and are unknowns that need to be 
l q l 

However, the retardation 


time. 

determined from experimental data. 
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time, X^ can be forced to be the same for each I th Kelvin 
element in each of the four material directions. This is 
reasonable since is predetermined or fixed when fitting a 
Prony series to experimental data with only allowed to 

vary. 


Furthermore, each layer in a composite laminate will 

have a set of four S terms describing its compliance 

<1 

matrix. If all the layers are of the same type of material 

and not rotated, i . e. all 0° direction, then ,S * S where 

k q q 

k is the ply layer in the laminate and Eq. S. 3 becomes 

Li s , ,L i s * * 
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However, if ply k is rotated, then the compliance matrix 
becomes fully populated, 
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Where each of the , ,S .. terms can be calculated from the 

kl qtj 

transformation matrices [19,20]. Similarly, the compliance 
coefficients, ,D , which will be used exclusively from this 
point on, can also be written in matrix form and rotated 
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where for fiber direction term, q = 1 ; 
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D ■ 4 m* nf ,,D 

kl aaa k k kl i 


for shear term, q = 4; 
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where m = cosCd D , n = sinC#}, D is the unrotated and 

Jc k k k kl q 

D . the rotated compliance terms C& is the angle of 

kl qij 

rotation}. Although it seems unnecessary and overly complex 

to split D .. into four parts, one for each material 
kl qtj 

property direction, this allows different stress nonlinear 
effects to be modeled in each of the four direction, which 
will be devel oped 1 ater . 

Unlike the compliance terms, ,D , the relaxation 

kl q 

times, X , are constrained to be the same in each of the 
four directions which eliminates the need to rotate them. 

All layers or plies are also assumed to be made of the same 
material which alleviates the need to kept track of the ply 
number when dealing with X . There are, however, some 
limitations on X . There should be at least one Kelvin 
element for every 1^ decades of time that is being examined 
since the effect of the Kelvin element is only felt over 
that time period. The common practice is one Kelvin 
element, thus one relaxation time, X , for every decade of 
time. For orthotropic materials it is further convenient to 
set X^ the same in all material property directions. A 
typical Prony series might have X^= 1, \^= 10, X ■ 100, 
etc. , for each ply and direction. 

If all the stresses in each layer and direction were 
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constant, then the time dependent strain could be easily 

calculated at this point by substituting the Prony series 

CEq. 5. 43 for each direction into the constitutive equations 

relating stress and strain CEq. 5.63 and solve for the 

desired time. However the stresses in each ply can in fact 

change with time which means stress is a function of the 

current strain rate as well as the current strain. Even 

though the restriction of constant stress was used to get 

the original Prony series in characterizing the material , 

that restriction is not necessary true in the actual 

numerical solution process. The matrix Eqs. 5.3, 5.5, and 

5.6 are still needed to show how the compliance terms can be 

manipulated and rotated to obtain the . ,D .. terms but they 

kl qtj 

are not used to obtain the strain. Instead the strain and 
stress equilibrium equation can be employed to calculate the 
strain. However a expression for the strain without the 
strain rate must first be developed. 

The original differential equation, Eq. 5.1, can be 
rewritten as 


k l 
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rt = < 
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kl i j 


kl^ij 
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where ,,e. . and . are the strain rate and strain, 

kl ij kl ij 

respectively, and , , a. . is the stress in each Kelvin 

kl ij 

element* l* ply* k* compliance direction* q* and rotated 
position* Ci*JD. This equation can be approximated by 


t+i t 

kl^ij kl^ij 
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where At is the time step size, t+1 is the new time and t is 
the old time. This particular approximation is called a 
Backward Euler Method CBEM3 and is classified as a first 
order implicit method. By using an implicit method, the 
solution, will be unconditionally stable regardless 

kl tj 
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of the time step size. This is not to say that it will 
converge to the correct answer but it will not diverge or 
blow up. This ’unconditionally stable’ characteristic only 
holds true for the first and second order implicit numerical 
approximations C6]. Higher order implicit methods and all 
explicit methods are only conditionally stable, i.e. has a 
maximum time step before it might diverge. 

The BEM, a first order implicit method, will be 
examined in detail in the remainder of this section. The 
second order implicit method, call the Backward Trapezoidal 
Method CBTMD is developed in appendix A. Equation 5.9 can 
be rearranged to give 
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where At = h. The unknowns are and cr t+i while h, 

kl ij kl ij 

X, , and ,,D .. are given and , ,£. . is known from the 

l kl qij kl ij 

previous time step. The total creep strain, for a 

particular direction, i, and layer, k, is simply the sum of 
all the creep strain in that direction 
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It should be noted that Eq. 5.10 is only for linear 
viscoelastic material. In order to include nonlinear stress 
effects, ,D needs to be modified to become a function of 

kl qtj 

stress. This is easily done by multiplying by a 

dimensionless stress function which would account for any 
nonlinear stress effects such as 
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where /CoO is a stress dependent function. This type of 

formulation allows the nonlinear compliance at any stress 

level to be scaled up or down from the linear compliance. ^ 

As an example, consider the linear compliance D represented 

by a single Kelvin element 

D* = D(l-e~ t/V ) C5.13D < 

If the nonlinear stress function is assumed to be / CoO = 

2 

Cl+aoO then the nonlinear compliance would be 

< 

d"= D(l-e _t ' /X ]*C:i+ao' 2 :> CS. 14Z) 

where a is a constant and a is the current stress. If the 

linear and nonlinear compliance curves are graphed, the ( 

2 

scaling factor, Cl+aoO, is quickly identified Csee Fig. 

5. ID. This formulation only works if the nonlinear stress 
can be described by a vertical shifting of the compliance 
curves. { 

Vertical shifting of the basic compliance curve to 
account for nonlinear stress effects is a common method of 
modeling nonlinear viscoelastic responseC 21 , 22] . Most 

nonlinear viscoelastic models such as the Schapery, Findley, 1 

and other power law based models employe this concept by 

using nonlinear stress functions. Figure 5.2 shows a simple 

nonlinear power law with a nonlinear stress function, /CoO , 

and how it is scalable. Since this vertical shifting 1 

concept works well for power law based models it should also 

work for a Prony series since, for many cases, the Prony 

series will just be a fitted equation to a Power law model. 

The Prony series can be scaled by just scaling the 1 

compliance coefficients, D, for a particular material 
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direction by the same amount that the power law would be 

scaled. Thus all Kelvin elements, l, will be have the same 

nonlinear scale factor, / C, oO . However, since the stresses 

q k 

are different in each ply, k, there will be a different 
scale factor for each ply. 

The nonlinear stress parameter, ^cr, used in the 
nonlinear stress function can be any function of the matrix 
or fiber stress states. A common parameter for the 
transverse and shear nonlinear compliance is the octahedral 
shear stress in the matrix which is a function of matrix 
transverse stress, a , and matrix shear stress, a A 

m2 m 12 

more detailed explanation of octahedral shear stress 
parameter can be found in the report by Dillard, et al Cl]. 
It is sufficient to say at this point that ^ a will be a 
function of the ply’s stress state, a , a , and or 

k 1 k 2 k 3 

regardless of the complexity. 

To introduce the nonlinear compliance function into the 
general formulation, substitute equation 5.12 into equation 
5. 10 to give 
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Where /CoO is evaluated at t+1. Note that /Co ,t * 1 } is a 

function of the future ply stress state, ,o ,t+1 , .o' 1 * 1 , and 

k l k 2 

.o' 1 * 1 , all of which are unknown. Thus Eq. 5.15 can be a 

k 9 

complex nonlinear function of a which necessitates the 
need for a numerical solution method. Equation 5. 15 can be 
rewritten as 
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Where uSti - CX^+hT kl D qij 

„t _ kl^ij *"1 
kl ij C^+h} 

t+i t+1 

a — a 
k j kl ij 

The subscript, l has been dropped from a because the stress 
is the same for all Kelvin elements in a particular material 
direction since they are all in series. The subscript i was 
also dropped since stress is a vector and not a tensor 
quantity. Similarly, the subscript j was dropped from the 
strain. The dropping of subscripts i and j for the stress 
and strain quantities, respectively, can be understood by 
reviewing the matrix equations Eq. 5.2 and 5.3. 

Substituting equation 5.15 into 5.11 will give the 
total creep strain for ply, k , and direction, i. 
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For the total strain, the elastic strain also needs to be 
added to equation 5.14. The elastic strain can be modeled 
as a nonlinear spring in series with the Prony series for 
each direction. 
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q=l 


where D is the linear compliance of the spring Cl=0}, 
ko qij 
t + l 

g C a } is the nonlinear function of stress for each 

q k « t*i 

direction Csimilar to /Co-}}, and , e. is the elastic 

ko v 

nonlinear strain. Adding equation 5.17 and 5.18 gives 


Page 30 


t t+i 
k*\ 


it ■ f L c ,ij 

j = l 1=1 q= 1 ^ J 


E 

k l i j 


I r*N ^ t ♦ 1^ I t+i 

) L d # c u a 3 

/ . IkO qi j q k Ik j 

q=l ^ J 


or 


—-it 

j = 1 1 = 1 q = 1 L J 

* qCk ° ,w ' : ’ -<1 M*] 


C5. 19} 


T tf 1 T t+i 

where c = c since it is assumed that layers deform 

k i v 

equally without debonding or damage. There are a total of 

T t+1 t+1 

3k+3 unknowns, c . and .o'. , but there are only 3k 

i. k j 

equations from Eq. 5.19. The additional 3 equations come 
from imposing stress equilibrium in each of the 3 stress 
directions 
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where °N is the input load on the laminate and N. is the 
J J 

actual load in each layer, k, when loaded. Equation 5.20 

can be rewritten in terms of stress to give 


o 



f t + 1 




C5. 21} 


where °a is the input stress and ,t is the thickness of 
j k 

each ply. This equation gives the 3 additional equations 
necessary to solve for the stress and strain unknowns at t+1 
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time. The equations are nonlinear in terms of stress and 
are ill conditioned. They can be solved by an iteration 
technique called the Newton-Raphson Method. Simpler direct 
methods such as Gauss -Seidel can’t be used since the 
coefficient matrix is not guaranteed to be diagonally 
positive. The Newton-Raphson takes longer to solve the 
nonlinear set of equations for each iteration, since the 
Jacobian matrix must be calculated, but it converges much 
more rapidly than the other direct iteration methods. 

In summary then, the three major problems with the 
current composite nonlinear viscoelastic analysis programs, 
nonlinear effects, stability, and orthotropic material, have 
been solved by using a differential equation formulation 
based on a series of nonlinear Kelvin elements. 
Implementation and results of this solutions technique are 
discussed in the next section. 
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6. Verification of the Nonlinear Differential 
Equation Method 


The method of solution presented in section 5, the 
nonlinear differential equation method CNDEMD, to calculate 
the nonlinear viscoelastic for orthotropic composite 
materials needs to be verified by comparing it to exact 
solution and other solution techniques. The section will 
present two simple examples, one based on the Kelvin element 
and the other on the power law, of a multilayered 
viscoelastic material for both linear and nonlinear cases. 
The solution will be compared to the exact solution, if 
obtainable, and other numerical solutions. 

The first example is a modified version of the example 
presented in section 3, that represents a simple 
one-dimensional two part material; one part is viscoelastic 
and the other elastic. The elastic material is modeled by a 
single spring and the viscoelastic material by a spring and 
a Kelvin element in series, as shown in Fig. 4.1. Since 
this example is relatively simply and one dimensional, it is 
possible to find a closed form solution for the linear case. 
For the nonlinear case, however, a Runge-Kutta method was 
employed to solve the resulting nonlinear first order 
differential equation. 

The linear case, assuming the applied stress, <y Q , is 
constant has a closed form solution of 
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The numerical results of the NDEM technique, both by the 
Backward Euler Method CBEMD and Backward Trapezoidal Method 
CBTMD, are shown in Fig. 6.2 along with the exact solution 
and the VI SLAP program technique solution. The spring and 
dashpot constants were assumed to be E^ E 2 = 1, E g = 0.11, 
tu = 1, and <y Q is constant, equal to 1, for all solution 
methods. The second order BTM solution matches the exact 
results closely whereas the first order VI SLAP and BEM 
solutions are high and low, respectively. This deviation 
can be accounted for by being only a first order solution 
technique. It is interesting to point out that the VI SLAP 
solution begins to oscillate and become unstable, as would 
be expected since it is an explicit solution method whereas 
the NDEM, for both BEM and BTM, is an implicit method. Also 
notice that the step size is large, 5 steps per decade, 
which would be considered the maximum step size but yet the 
second order NDEM is very accurate. 

The same basic model can be used for a nonlinear 
viscoelastic material by simply changing E g and y to include 
nonlinear stress effects. For the current nonlinear example 
E and y are as follows 

3 

a 

E = 0.1 — — + 0.1 C6. 2D 

a a 

2 

y * 10 | 0.1 + 0.1 J = 10 E g C6.3D 

where o^is the stress in material 2 CFig. 6. ID. This type 
of nonlinearity will cause the material to become stiffer as 
time progresses since the stress, o' , is decreasing in the 
nonlinear dashpot. As the stress decreases in the Kelvin 
element, the spring becomes stiffer and can ultimately carry 
more of the total load. Likewise, the nonlinear dashpot 
will become more viscous and the viscoelastic response will 
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be retarded. The other parameters are similar to the linear 
case, E = E = 1 and at = 1 . The results of both the NDEM and 

12 O 

VISLAP techniques are shown in Fig. 6.3 for the nonlinear 
Kelvin element. To obtain the exact solution one must solve 
a nonlinear equation of the form 

£ + K£ + Kc*+K=0 C 6. 4D 

12 3 

with tfCCD = c 

This equation is difficult to solve for a closed form 
solution but good results can be obtained by using a 
Runge-Kutta numerical method with small time steps. The 
results from a Runge-Kutta solution is plotted on Fig. 6.3. 
For a simple nonlinear example model, like the one being 
examined, it is possible to use a Runge-Kutta solution as a 
check which is a well proven and reliable numerical method. 
For the more general orthotropic problems such a method is 
not possible as discussed in the proceeding sections. 

Similar to the linear case, the first order solution 
methods, VISLAP and NDEM using BEM, are not as accurate as 
the second order NDEM using BTM. Also the explicit method, 
VISLAP, becomes unstable at long time steps. 

The second example case is again a two part 
one-dimensional material with one part viscoelastic and the 
other elastic. This viscoelastic material is modeled as a 
power law and a spring in series, and the elastic material 
as a spring. Figure 6. 4 shows the mechanical model 
describing this test case. 

The Power law parameters used are m = O. 1 and n = 0.25 

for the linear case and 

ct 

m = 0.1 — — + 0.1 C6. 5D 

a 

2 

and n = 0.25 for the nonlinear case. For both the linear 
and nonlinear case E^= E^= 1-0. The results comparing just 
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the VI SLAP and NDEM CBEM and BTMD methods are shown in Figs. 
6.5 and 6.6. Both the linear and nonlinear cases show the 
results for all methods very close, with the BTM between the 
other two methods. This is similar to the results of the 
linear and nonlinear Kelvin cases discussed earlier. There 
is no exact solution available to compare results and the 
resulting equation can’t be solved by the Runge-Kutta in a 
convenient manner. However the results of both the VI SLAP 
and NDEM techniques are similar, givinng some reassurance 
that the answer is correct. 

In summary, the nonlinear differential equation method 
CNDEMD in solving nonlinear viscoelastic problems that 
involve multiple material layers has been shown to be an 
accurate method and does converge to the correct answer. 

The two test cases examined, Kelvin and Power law models, 
showed the NDEM results match the exact solution and/or 
other numerical methods. The second order BTM technique 
proved to be the most accurate and was stable for all time 
steps . 
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7. Summary and Conclusions 

This report has looked at various methods to solve 
nonlinear viscoelastic problems that deal with orthotropic 
materials such as fiber reinforced composites. Earlier 
methods, such as the VISLAP computer program algorithm, was 
examined and some of the deficiencies discussed. The main 
three problems of these methods were ID, stability of the 
solution technique, 2D, time step size stability, and 3D, 
solution time length and computer memory storage. Two other 
methods were examined in detail , Volt err a Integral and the 
Zienkiewicz method, plus a new method, the Nonlinear 
Differential Equation Method CNDEMD was developed to try to 
over come some of the deficiencies. 

The Volt err a Integral allowed the implementation or 
higher order solution techniques but it had difficulties on 
solving singular and weakly singular compliance functions. 
The power law compliance function, which is weakly singular, 
was solvable only with very small time steps. This method 
also needs an every increasing amount of computer time as 
the solution process goes further out in time, similar to 
the VISLAP method. This was due to the hereditary type 
integral solution process which must recalculate the total 
integral for each addition time step. This method was found 
to be unacceptable for reasons of computer time needed and 
accuracy. 

The second method examined was the Zeinkeiwicz solution 
technique which requires the viscoelastic response to be 
modeled by a Prony series. This method works well for 
linear viscoelastic isotropic materials and small time 
steps. The biggest advantage of this technique is that the 
solution algorithm can be written in a recursive fashion 
which does not require the recalculation of the past results 
like the VISLAP and Volterra Integral methods. This allows 
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the solution at long times to be done efficiently and 
quickly. The biggest problem with this method is the limit 
on time step size since the method uses an explicit solution 
technique. Thus the solution can become unstable and 
diverge from the correct answer. 

To overcome the above deficiencies a new method* NDEM* 
was developed. This method requires the viscoelastic 
response be described by a modified Prony series which 
allows nonlinear stress effects to be included. The 
differential equations that model each of the Kelvin 
elements in the Prony series* are then solved 
simultaneously. By using the basic differential equations* 
an implicit solution method can be used. This causes the 
solution process to be unconditionally stable for any time 
step. The general method of solving the nonlinear 
simultaneous equation used was the Newton -Raphson method 
which assures convergence even if the coefficient matrix of 
the equations is not positive definite. In addition to 
overcoming the numerical problems this method was extended 
to include orthotropic nonlinear viscoelastic materials. 

The NDEM technique was shown to be accurate and stable 
on two test cases* Kelvin and Power law based* for both 
linear and nonlinear conditions. The advantages of NDEM is 
that it is stable for all time step sizes* the solution 
algorithm is stable and converges to the correct solution, 
and the computer time is minimized. 
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Appendix A - Backward Trapezoidal Method 


In Chapter 5 the Backward Euler Method was used in 
solving the nonlinear viscoelastic problem of orthotropic 
composite laminates and a detail derivation was given. The 
Backward Trapezoidal Method CBTMD will be briefly developed 
in this appendi x . 

Recall the basic differential equation, Eq. 5.8, of a 
single Kelvin element 

n'ii “ Z 

q= i 

where . £. . and . ^ are the strain rate and strain, 
kl ij kl ij 

respectively, , a. is the stress in each Kelvin element, l, 
k l tj 

ply, k, compliance direction, q, and rotated position, 
Ci,JD. Using the BTM the numerical approximation becomes 


D 

kl qtj 


a 

kl i j 


£ 

kl Ij 


CA. ID 



where At is the time step size, t+1 is the new time and t i 
the old time. This can be rewritten as 


t + 1 
kl^L j 


2\, +h _ 
l q= 1 


f t+i h r* t 

2 k l qij kl°ij 2\ +h 2 kl qij kl°ij 

— 1 l q = 1 


2\, -h 

l t 

2\^+h k l C i j 


CA. 3D 


where h = At. When the nonlinear stress function are 


< 
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included CEq 5. 103 Eq A. 3 becomes 

1 + 4 h t d + 

Z k l qij kl ij y q k 


£ 

kl i j 


2\ +h _ 
l q= 1 


ax -h 


ax +h _ 

l q=l 


7 l1 D . . a 1 . / C V) + pi l . K . .tf 1 . CA. 43 

2L,kl qij kl tj ; q k 2\ +h k l i J 


This can be further simplified as 


t+i 

ki C i 


yr c / cv +i d] v 

^[kl qij q k Jk j 

q=l 


1+1 + E l 
k l i j 
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where 


kl^ " 


y\c / cV) ] v 

/ . [_kl qtj qk Jkj 


+ 2 V h t 

2\ +h k l ^i j 


q= i 

h 


D 


kl qij C2X^+hD kl qij 


t+l t+1 

a = a 
k j kl ij 


Summing all the creep strains together for each element in 
the series of Kelvin element will give, similar to BEM, 


a 

n 

4 


•< 

c t +1 ^ 

«.'! ~L 

T 

I 

_ , . t + 1 ^ 

. . c . . / C a 3 

k l qt j q k 

t+1 . t 

Cr + E 
k j kl i j 

J = 1 

1 = 1 

\o | 
n j 

* - 
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From this point the derivation is the same as the BEM. As 
with the BEM, the BTM is unconditionally stable for all time 
steps. It is a second order method which will be more 
accurate than the BEM. 


i 
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Figure 3.2 Numerical Solutions of a 3 Parameter 
Dashpot Model Using the Volterra Integral 
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Figure 3.3 Numerical Solutions of a 3 Parameter 
Power Law Model Using the Volterra Integral 
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Figure 3.4 Numerical Solutions of the Volterra Integral 
Using the Simpson Rule with the 3/8 Rule 
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Figure 5.1 Vertical Scaling of a Single Kelvin Element 
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Figure 5.2 Vertical Scaling of a Power Law Compliance Function 
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Figure 6.3 Results of Various Numerical Techniques 
for the Nonlinear Kelvin Element Test Case 
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Figure 6.5 Results of Various Numerical Techniques 
for the Linear Power Law Test Case 
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